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ABSTRACT 

Context. Molecular clouds typically consist of 3/4 H 2 , 1/4 He and traces of heavier elements. In an earlier work we showed that at 
very low temperatures and high densities, H 2 can be in a phase transition leading to the formation of ice clumps as large as comets or 
even planets. However, He has very different chemical properties and no phase transition is expected before H 2 in dense interstellar 
medium (ISM) conditions. The gravitational stability of fluid mixtures has been studied before, but these studies did not include a 
phase transition. 

Aims. We study the gravitational stability of binary fluid mixtures with special emphasis on when one component is in a phase 
transition. The numerical results are aimed at applications in molecular cloud conditions, but the theoretical results are more general. 
Methods. First, we study the gravitational stability of van der Waals fluid mixtures using linearized analysis and examine virial equi¬ 
librium conditions using the Lennard-Jones intermolecular potential. Then, combining the Lennard-Jones and gravitational potentials, 
the non-linear dynamics of fluid mixtures are studied via computer simulations using the molecular dynamics code LAMMPS. 
Results. Along with the classical, ideal-gas Jeans instability criterion, a fluid mixture is always gravitationally unstable if it is in 
a phase transition because compression does not increase pressure. However, the condensed phase fraction increases. In unstable 
situations the species can separate: in some conditions He precipitates faster than H 2 , while in other conditions the converse occurs. 
Also, for an initial gas phase collapse the geometry is essential. Contrary to spherical or filamentary collapses, sheet-like collapses 
starting below 15 K easily reach H 2 condensation conditions because then they are fastest and both the increase of heating and opacity 
are limited. 

Conclusions. Depending on density, temperature and mass, either rocky H 2 planetoids, or gaseous He planetoids form. H 2 planetoids 
are favoured by high density, low temperature and low mass, while He planetoids need more mass and can form at temperature well 
above the critical value. 

Key words. Instabilities - ISM: clouds - ISM: kinematics and dynamics - ISM: molecules - Methods: analytical - Methods: 
numerical 


1. Introduction 

Typically, the Milky Way molecular clouds consist of molecular 
hydrogen (^H 2 ) and helium (^He) in the respective mass fraction 
of ~74% and ~24% and traces of heavier elements in the form of 
atoms, molecules, and dust grains (Draine 2011) . The He mass 
fraction is thus non-negligible. Even though H 2 and He are by 
far the most abundant chemical components, they remain hardly 
detectable, and most of the time they are inferred from CO emis¬ 
sions (Bolatto et al. 2013). Thus, the dynamical and chemical 
processes associated with H 2 and He in molecular clouds are 
still poorly known, especially when considering sub-AU scales. 

In Fuglistaler & Pfenniger (2015, hereafter FP2015), we 
discussed substellar fragmentation including gravity in single 
species fluids presenting a phase transition, such as very cold 
molecular hydrogen in molecular cloud conditions. We showed 
that fluids in a phase transition (i.e. subject to a chemical in¬ 
stability) are anyway also gravitationally unstable because any 
density fluctuation is not compensated by a pressure variation, 
but by a change in condensed matter fraction. In phase transi¬ 
tion conditions arbitrary small condensed clumps can form. The 
possibility of forming H 2 ice clumps in the ISM, from grains, 
to comet-like bodies to rocky or gaseous planet-like bodies pro¬ 


vides a scenario for baryonic dark matter extending the scenario 
of Pfenniger et al. (1994); Pfenniger & Combes (1994) towards 
micro-AU scales. However, since molecular clouds contain a 
substantial fraction of He, it is necessary to investigate how this 
component might modify the findings of our previous study. 

Although at first sight from the chemical point of view both 
H 2 and He present an outer electronic shell made of two elec¬ 
trons, their chemical properties differ markedly, mainly because 
of quantum physics. The individual properties of H 2 and He are 
well known from laboratory data (Air Liquide 1976) and shown 
in Fig. 1. H 2 and He are in a phase transition when on the con¬ 
densation wall linking the gaseous and solid or liquid phase. He 
has a lower critical temperature than H 2 (5.2 K vs. 32.9 K) and a 
lower critical pressure (0.227 MPa vs. 1.286 MPa). In the highly 
dynamical conditions present in molecular clouds, such as su¬ 
personic turbulence (Elmegreen & Scalo 2004) , phase transition 
conditions may be reached thanks to a combination of pressure 
increase and/or temperature decrease. In such a case, phase tran¬ 
sition conditions are reached for H 2 well before He. The con¬ 
ditions of phase transition of the mixture H 2 -He may, however, 
change the conclusions made in the single species case. 

The properties of H 2 -He mixtures has been mostly studied 
in detail for temperatures above the critical and at high densi- 
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ties (e.g. Streett 1973; Koci et al. 2007; Becker et al. 2014) and 
especially in conditions similar to gas giant planets (Vorberger 
et al. 2007; Saumon et al. 1995) . Taking quantum effects into 
account, Safa & Pfenniger (2008) calculate the thermodynamic 
properties of H 2 -He mixtures below critical temperature from 
the known intermolecular potentials and obtain the critical point 
and stability of the mixture itself. 

The gravitational stability of self-gravitating binary and mul¬ 
ticomponent fluids has been studied by Grishchuk & Zeldovich 
(1981) , who showed that there can be only one unstable solution. 
If a fluid mixture is gravitationally unstable then all components 
are affected. They note that in the case (dPIdp)^ = 0 (i.e. the 
sound-velocity formally vanishes), which is the case in a phase 
transition, the fluid is always gravitationally unstable, but they 
do not go into more detail on that specific case. Jog & Solomon 
(1984a,b) first discussed the stability of two-component disks. In 
a similar fashion de Carvalho & Macedo (1995) studied oscilla¬ 
tions and resonances in a binary fluid mixture. Volkov & Ortega 
(2000) discussed the stability of self-gravitating systems with a 
spectrum of particle masses and consider rotating mediums. 

When a phase transition occurs in the presence of external or 
internal gravity, the fluid dense phase may precipitate in the form 
of rain, snow, or hail in the atmosphere, leading to a fragmen¬ 
tation that is impossible to describe with usual hydrodynamic 
codes, in which a single phase in local thermal equilibrium is 
implicitly assumed. In FP2015 we showed that method phase 
transition and precipitation can be simulated for a single species 
with molecular dynamics. The possible objects condensing from 
the gaseous phase can take various masses, typically covering 
the entire range from grains, comets to planets or larger. With 
two species with different molecular weights the number of pre¬ 
cipitation scenarios that can be envisioned increases. Could it 
be, for example, that bodies form with a core made of solid H 2 
surrounded with an atmosphere of H 2 and He or that a solid H 2 
crust floats on a gaseous He core? 

To answer such questions we use the same molecular dynam¬ 
ics code as in FP2015 just adding a second species, and scaling 
the particles properties to the respective properties of H 2 and 
He. We control the finite number resolution effects by perform¬ 
ing simulations over a range from 1.25 • 10^ to 80 • 10^ particles. 
We restrict the investigations to the simplest set-up combining 
gravity with molecular dynamics. To control gravitational in¬ 
stability, we investigate a single plane-parallel collapse in one 
direction of a periodic cube, where the initial temperature and 
density are simulation parameters. As explained in Sect. 2.2 and 
Appendix B, the collapse geometry (sheet-, filament-, or point¬ 
like) is crucial to reach phase transition conditions starting from 
typical ISM conditions. Sheet-like collapses (pancakes) can in¬ 
deed lead temporarily to very dense conditions without much 
heating, contrary to the other cases. 



Fig. 1: H 2 (blue) and He (red) phase diagrams. For clarity, only 
a part of the upper, almost constant density condensed phases of 
both species and the low-density gas phase of He are shown. 


to the residual kinetic energy when the bulk translational, ex- 
pansional, and rotational velocities are subtracted, be it glob¬ 
ally or locally. Strictly, this definition is operational and useful 
only if the velocity distribution is unimodal and its second or¬ 
der moment exists. Further detailed discussion about this topic 
would be out of scope, as in this article we consider either global 
or local temperatures for particle systems with no or negligible 
amount of ordered motion, so the stated temperature is equiva¬ 
lent to the particle kinetic energy. The high degree of collisional- 
ity in molecular interactions ensures the rapid destruction of any 
initial correlations leading to the convergence towards thermal 
states. 


2.1. Jeans instability 

Considering a fluid mixture as a one-component fluid using av¬ 
erage quantities such as the density p and pressure P, the Jeans 
instability criterion (see App. A.l) would be 


,2 < , 2 ^ 

^ {dPjdpX 


( 1 ) 


2. Gravitational stability of a fluid mixture 

In a fluid consisting of K different components /, the total mix¬ 
ture number density is ^ the mixture mass density is 

P - Y^iPi with pi = miHi, and the mixture pressure is P = Pi- 
In the case of an ideal gas, Dalton’s law states P = 2/ ^iPi- Each 
component has a molecular fraction xt = niln and a mass frac¬ 
tion Wi = nii/m with m = 

The notion of global temperature in a system with long-range 
forces is an unsettled topic as the key assumption of extensiv- 
ity in thermodynamics breaks down in long-range force systems 
(e.g. Padmanabhan 1990 ). When dealing with particle systems 
we can however always define the temperature as proportional 


where k is the wavenumber, k] the critical Jeans wavenumber, 
and G the gravitational constant. This is, however, inappropriate 
if each component has a different mean square velocity, which is 
the case for an isothermal fluid with different molecular masses. 

In order to correctly predict the stability of a many- 
component fluid mixture with different densities pi and partial 
pressure P/, each component has to be treated individually: the 
Grishchuk-Zeldovich criterion (see App. A.2), which is the sum 
of each component Jeans’ criteria, reads 


k^ < k 


2 ^ 

GZ - 


Z 47rGpi 
dPildpi ■ 


( 2 ) 


Article number, page 2 of 22 





















A. Fuglistaler and D. Pfenniger: Formation of H 2 -F[e substellar bodies in cold conditions: 


2.1.1. Ideal gas mixture 


A fluid far from the condensed phase can be approximated with 
the ideal gas law 


p 

= nk^T, 

dP\ 

ksT 


= 7 - 

m 


(3) 

(4) 


where is the Boltzmann constant, y the adiabatic index, and 
m the molecular mass. The partial pressures can be calculated 
using Dalton’s law. Equ. (2) becomes 


- 

^GZ,id “ 


47rG 

yk^T 




(5) 


where id stands for ideal gas. This equation differs from Equ. 
(1), which in the ideal gas case becomes 


^J,ic 


47rG 

ykBT 


nm 


( 6 ) 


In the case where all components have the same temperature 
fcz,id ^ ^J,id independent of the molecular fractions Xi = nijn 
and molecular mass m^. In the case of a H 2 -He mixture, the max¬ 
imum k^ziJ^nd ~ ■^H 2 = 0.61. 


2.1.2. van der Waals fluid mixture 


When approaching a phase transition, the ideal gas law does not 
take into account condensations and does not yield correct values 
anymore. We showed in EP2015 that the van der Waals equation 
of state (van der Waals 1910) describes a phase transition rather 
well provided that the Maxwell construct is taken into account 
(Clerk-Maxwell 1875; Johnston 2014), i.e. 


_ / 24r, 6nr\ 

\ dp ^\m{nr - 3)2 m / ’ 


(7) 

( 8 ) 


in gaseous and solid/liquid form with the reduced values = 
PI Pc, Tr = r/Tc, = nj He and the critical pressure, temperature 
and density P^ Tc and Uc. In the case of a phase transition, P^ = 
const (Maxwell construct) and (dPrldp)s = 0. 

The Maxwell line is very similar to the laboratory condensa¬ 
tion line for H 2 in a T - P diagram as can be seen in Pig. but is 
rather off for He, especially at low temperatures. As in the astro- 
physical context, correctly representing the H 2 phase transition 
is essential for our study. A H 2 phase transition always occurs at 
a lower pressure-temperature ratio than for He. 

In the phase transition regime (dPldp)s = 0, varying density 
allows the pressure to remain constant. This is a crucial prop¬ 
erty for this work, since gravitational contraction is no longer 
compensated by pressure increase. We show in App. A.2 that a 
two-component fluid is always gravitationally unstable as soon 
as (dPildpi)s = 0 for any component i. In a similar fashion, the 
same can be deduced for ^-component fluids (Grishchuk & Zel- 
dovich 1981). 



Pig. 2: H 2 and He laboratory data and van der Waals vapour 
curve derived from Maxwell construct. Adiabatic compression 
curves of an initial sphere to sheet-, filament- and point-like ge¬ 
ometries for interstellar conditions (T = 10 K, P = 10“^^ Pa) are 
shown, as explained in Appendix B. Sheet-like collapses are al¬ 
lowed to reach the phase transition regime even without cooling. 
Cooling would displaces the curves to lower temperature. 


In addition, as shown in App. B.l, the adiabatic matter com¬ 
pression of a sheet-like geometry leads to a finite increase of 
potential energy, leading to a maximum relative temperature in¬ 
crease of only 2.1, while the energy diverges logarithmically 
in a filament-like geometry and as in a point-like geom¬ 
etry, where Z = Pfinai/Pinmai is the density compression. This 
can be seen in Pig. 2^ which shows how a sphere moves in this 
diagram when adiabatically compressed towards a sheet, fila¬ 
ment, or point initially in interstellar conditions (T = 10 K, 
P = 10“^^ Pa). Whereas the temperature is quickly increasing 
with filament- and point-like geometries, in the sheet-like geom¬ 
etry it only rises to ~21 K, which is well below the 33 K critical 
temperature of H 2 . 

When taking radiative cooling into account, the temperature 
increase by contraction is even smaller. In App. B.2 we show 
that the opacity of a sheet-like collapse is barely increasing. If 
the initial medium is transparent, the final sheet is also transpar¬ 
ent. On the other hand, in filament- and point-like collapses the 
opacity increases approximately as a power 1/2 or 2/3 of com¬ 
pression, quickly reaching a full opacity regime able to stop the 
collapse. 


2.3. Lennard-Jones mixtures 


2.2. Plane-parallel collapse 

Pin et al. (19^) and ZeTdovich (1970) show that a plane- 
parallel collapse, leading to a sheet-like geometry, is a faster col¬ 
lapse than filament- or point-like geometries. This has been con¬ 
firmed using numerical simulations by Shandarin et al. (1995). 


The Jeans instability of Equ. (2) requires thermal equilibrium 
and is simplified by only considering linear perturbations. It does 
not predict its non-linear evolution when unstable. This is a mo¬ 
tive to use molecular dynamical simulations with a Pennard- 
Jones potential in addition to gravity for studying such non¬ 
linear phenomena. 
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The Lennard-Jones potential 


tuCr) = 4- 
m 



(9) 


reproduces the van der Waals equation of state when in equilib¬ 
rium conditions using the following relations (Caillol 1998): 


Tc 

= 1.326-^, 

(10) 

tic 

= 0.316cr-V 

(11) 

Pc 


(12) 


For a fluid mixture, we use the usual Lorentz-Berthelot combin¬ 
ing rule for the Lennard-Jones potential between two molecules 


'I>Lj(r,7) = 4^ 

rui 



(13) 


where (Tij = ^(cri-\-crj) and eij = (Lorentz 1881; Berthelot 
1898) . The most accurate mixing rules tor energy and distance 
are (Banaszak et al. 1995; Chen et al. 2001) 


K K 


i i 

(14) 

ZfZfxiXjecrl 

€ = -z- . 

cr^ 

(15) 


Since oc e and ric oc cr~^, we have ^c/^c,Qr = (o'/o'a)~^ and 
TclTc,a = Using Equ. (7), we find Pc = (8/3)rc^c, and 
therefore PJPc.a = (crlcra)~^ ■{el£a)- 

As in the ideal gas case, using these mixing properties to cal¬ 
culate k] in analogy to Equ. (1) instead of using Equ. (2) would 
lead to different results. In the ideal gas case, the biggest differ¬ 
ence is -10%, but in the present case of a van der Waals fluid the 
difference can be much bigger, for example k] always returns a 
finite number if the critical temperature of the mixture Tc > 1, 
whereas ^gz can nevertheless be infinite if one of the compo¬ 
nents is in a phase transition. 


2.3.1. Binary mixture 

Considering a fluid consisting of two components a and /J, we 
define 


TcJ3 _ 


(16) 

Pc,a 

/ \ 

_ 

- ’ 

(17) 

^c,a 



mp 

rua ■ 


(18) 


The following Lennard-Jones properties can be derived using 
Equ. (23-28): 



Eig. 3: van der Waals phase diagram with Maxwell construct for 
the H 2 - He mixture, each species considered as independent. 


and the following van der Waals mixture properties: 


^C,X 

^c,a 

= 


^a^P / 

■ 4 1 

Tc,^ 

Tc,a 

= 

^c,x 

^c,a 


Pc,X 


^c,x 

T'c.X 

Pc,a 


^c,a 

Tc,a ’ 


' y 




gm 


/ ^ ^ xi0 


( 21 ) 

( 22 ) 

(23) 


where xp = 1 - Xa. Without loss of generality we choose 6 < 1. 

Three different cases can be distinguished, as shown in Eig. 
3: 


(A) Tc^js < Tc,a < T 

Having the temperature above both critical values, neither 
component can be in a phase transition and kcz is always 
finite. 

(B) Tcfi < T < Tc^a 

P is still above the critical temperature and cannot be in a 
phase transition, but a may or may not be in a phase transi¬ 
tion depending on Ua. 

(C) T < Tc^p < Tc,a 

Having the temperature below both critical temperatures, 
both components can be in a phase transition. At equal com¬ 
ponent number density, a is faster in a phase transition, but 
theoretically, at low np, (3 could be in a phase transition with¬ 
out a, but this hardly ever happens in reality. 

These three cases are studied using computer simulations (see 
Table 1). 


O-aj} = ^Cr„(l + V , 

^ap = , 


2.3.2. Hydrogen-helium mixture 

(19) 

The critical temperature of H 2 and He are 32.97 and 5.19K, 

(20) whereas their usual Lennard-Jones e values are 36.4 and 
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ture. 


10.57 ^bK. This is in conflict with the temperature conversion of 
Equ. (16), as ^H 2 -He = 6.35 using critical temperatures whereas 
^H 2 -He”= 3.44 using the Lennard-Jones e values. The same is the 
case to a lesser degree for the critical density. 

All performed simulations are molecule independent, but 6, 
y, and need to be deflned. These values were set using Tc, 
^c, and the molecular mass of laboratory He and H 2 data (Air 
Liquide 1976) . The goal of this article is to understand the role 
of a secondary component in a fluid presenting a phase transition 
together with gravity. In molecular clouds, the most likely case 
of a phase transition is > T > Tc,uq, thus it is important to 
have a correct (Tc,He/^c,H 2 ) fraction. 





Fig. 5: van der Waals phase transition and unvirializable density 
for a binary mixture with Xa = 0.75. 

phase. Figure £ shows R and A as functions of an abundance 
number fraction. 

The above terms are for a fluid with no spatial separation 
of the species, i.e. a fluid in its initial state. The terms predict 
how an unstable fluid evolves. However, once a phase transition 
or a gravitational collapse happens, the species may separate. 
In that case, the terms have to be calculated for each species 
independently, using = I and 0, respectively. 


2.3.3. Virial theorem 

In FP2015, the Fennard-Jones potential has been decomposed 
in attractive and repulsive terms: Olj = Oa + In a binary 
mixture, we have to further distinguish the one-component terms 
Oq ,2 and Oyj 2 , and the cross-terms and (see Equ. 13). 

In analogy with Equ. (5) of FP2015, the virial theorem be¬ 
comes as follows: 

0 = 2£^]Qn -|- 12£^i- -|- 6£^a F • (^4) 

>0 <0 


2.3.4. Unvirializable densities 

The sum of the Eennard-Jones and kinetic terms must be positive 
as the gravitational term is negative: £"0 < 0 < F’r + F’a + £kin- 
This is the case if 

2Aea {crlnf <k^T + ARe^ . ( 31 ) 

There can be no virial equilibrium for any mass M > 0, if the at¬ 
tractive Eennard-Jones term dominates the kinetic and repulsive 
terms. We deflne the unvirializable density domain as follows: 


There are two negative, attractive terms, and two positive, re¬ 
pulsive terms. If the attractive and repulsive terms equalize each 
other, the system is in virial equilibrium. 

In the case of a homogeneous density and species distribu¬ 
tion of mass M, the energy terms are 


^kin 

Er 

E, 

Eg 


2m 

R-^ —M, 
m 

EaCrtn^ 

-A ^ M, 
m 

-G /g (n, m) , 


(25) 

(26) 

(27) 

(28) 


with /g > 0 depending on the geometry. The repulsive and at¬ 
tractive constants are 


R 

A 


xl + ev-^xl + 


om 

212 

26 


xl + 0V + ^[l+V [xaXl + xlxp) 
(l + + xlxp) 


(29) 

(30) 


where the lattice coefficients and Ca depend on the specific 
crystal lattice (FP2015), which are of importance for the solid 


D = [n '. H- < n < n+} 

1± 


n, = cr 


4! 


(32) 


If the term in the square root is negative, there is no real solution 
and therefore no unvirializable densities. This is the case if 


T>Trr 


4RkB ■ 


(33) 


Figure _5 shows the n± values and the domain of the van der 
Waals phase transition for H 2 -He mixtures with a molecular frac¬ 
tion of Vq, = 0.75. There are unvirializable densities above the 
critical temperature up to r^ax even though no phase transi¬ 
tion is possible. 

The domain of phase transition does not change a lot for the 
different xu^ > 0, as the H 2 phase transition temperature remains 
the same and the density changes as = v • ^. It is at much 
lower temperatures for xu 2 = 0, as there is no H 2 anymore and 
the phase transition domain switches to He. 

The evolution of fluids with unvirializable densities depends 
whether they are in a phase transition or not. If these fluids are in 
a phase transition, there is a gravitational instability independent 
of M (see 2.1.2), which leads to a collapse and the formation of 
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bodies of small mass. If they are not in a phase transition, there 
is only a gravitational collapse above a certain mass M (see Equ. 
^). Below that mass, clumps form, which leads to an augmenta¬ 
tion of the kinetic and repulsive Lennard-Jones terms until Equ. 
(31) is fulfilled, at which point an equilibrium is reached and the 
fluid may remain stable. 


first, as explained in Sect. 2.2, a velocity perturbation in form of 
a small plane sinusoidal wave in the v direction is superposed 
to the fiuid’s Maxwellian velocity distribution. The perturbation 
strength is of 1%; see EP2015 for how the perturbation is calcu¬ 
lated. 


2.3.5. Dynamical friction 

When discussing gravitational collapses, the concept of dynam¬ 
ical friction (Chandrasekhar 1943) is important, since heavy ob¬ 
jects may condensate from the gas and start to precipitate, i.e. 
move with respect to the gas in the local gravity field. Consider¬ 
ing a uniform density and Maxwellian velocity distribution, the 
dynamical friction of a heavy object of mass mh reads 


3.1. Units 

All simulations are performed in dimensionless units and only 
the ratios of physical quantities matter. The initial properties of 
a fiuid are the ratios of its temperature to the critical value Tc, its 
number density to the critical value Uc, the Lennard-Jones con¬ 
stant ratios 6, v, the mass ratio p, the molecular fraction Xq, and 
the gravitational potential strength yq. The needed molecule pa¬ 
rameters were set in accordance to laboratory H 2 and He values 


dvh 

dt 


47rG^mh p log(A) 


erf (X) - 


2X 


exp (-X^) 




(34) 


where log(A) is the Coulomb logarithm and X = v\^l V2cr„. For 
the qualitative analysis needed in this work, it is enough to state 


dt 


oc -mh . 


(35) 


This can be used in different cases. Eirst, He is twice as heavy as 
H 2 and can therefore be considered a heavy object and in some 
conditions lead to sediment faster than H 2 . Secondly, when H 2 is 
in a phase transition, H 2 ice grains may sediment faster than He. 


3. Method 


=0.157, 

(41) 

^C,H 2 


77 c He 

’ = 1 . 12 , 

(42) 

^C,H 2 


’^=2. 

(43) 

mu2 



The time unit is defined as the particle crossing time for the box 
of length L, i.e. 



(44) 


where = ^Nk^T j 2/ i=l...N. 

The gravitational constant strength is measured by a factor 
7 j relative to the ideal gas Jeans limit strength Gj, defined as 


Eor all of the simulations, the Large-Scale Atomic/Molecular 
Massively Parallel Simulator (LAMMPS) is used (Plimpton 
1995). The use of its short-range Lennard-Jones solver and long- 
range gravitational solver, super-molecules, and the rRESPA 
time integrator are discussed in EP2015. 

We recall the following super-molecule properties, where 77 
is the number of molecules per super-molecule. 


^SM 

= 77777 , 

(36) 

^SM 

= ye, 

(37) 

O' SM 

II 

(38) 


where m, cr, and e are the values for one molecule. The gravita¬ 
tional constant, described in molecular dynamics units (cr = e = 
m = 1), is 


Gsm 


Gm^ 

ere 


^2/3 


(39) 


In order for the gravitational force between two super-molecules 
to be consistently small compared to the intermolecular forces, 
77 should satisfy the following constraint: 


2 

773 


/ _5 

Gm2 V ^ 



(40) 


where is the cut-off radius in cr units (set to Tc = 4 in the 
simulations). 

Two molecules are considered as bound in LAMMPS if their 
distance is smaller than 1.3625(T. A clump of bound molecules 
can be either gaseous or condensed, dependent whether its tem¬ 
perature is above or below the critical value. 

Initially, the fiuid is uniformly distributed in a periodic cu¬ 
bic box. To reproduce the most generic plane-parallel collapse 


n 


Gj 


_G 

TrykBT 
2 njrrP: 

j=a,p ^ 


(45) 


3.2. Visualization 

In order to visualize the particle snapshots, two-dimensional 
number density maps are used. The introduced perturbation is 
along the x-axis, thus the sheet-like collapse is parallel to the yz- 
plane. Lor that reason, the density map shows the y- and z-axes 
where most of the relevant events can be observed. When show¬ 
ing all N particles, smaller aggregates are washed out and barely 
visible, which is why only a slice of the whole simulation box 
is shown. Ligure 6 shows how this slice is selected: The highest 
number density is determined in the x direction in order to be 
centred around the collapse region. Prom there, the slice width 
Ax is calculated in order to contain AsUce = /slice A particles. In 
that way, the slice width Ax differs for every snapshot, but always 
contains the same number AsUce of particles. 

The number density n and number fraction Xq, are represented 
by brightness and colour. As at low density the brightness is 
maximum and the colour is simply white, the colour map is best 
visualized in polar coordinates, where n is the radius and Xq, the 
angle. Ligure 7 shows the colour map used, for better contrasts, 
the brightness is represented in logarithmic scale. 


4. Simulations 

In PP2015, we simulated fluids with only one component. We in¬ 
troduced the terms comets for clumps that are principally bound 
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Table 1: Simulation parameters 


Name 



n 

T 


71 

Mot 


L/Re 

a b 

r/kyrs"' 



[«c,a] 

[m-2]« 

[Tc,a\ 

[KY 



(f) 

(g) 

(f) 

(g) 

(f) 

(g) 

AlO 

1.0 







0.34 

1.7 

27 

47 

0.86 

1.5 

A75 

0.75 







0.094 

0.49 

16 

28 

0.58 

1.0 

A50 

0.5 

10-2 

9.3 • 1025 

1.5 

49.5 

0, 0.5, 1.5 

1005 

0.039 

0.2 

11 

20 

0.45 

0.78 

A25 

0.25 







0.019 

0.1 

8.6 

15 

0.36 

0.64 

AGO 

0.0 







0.011 

0.056 

6.8 

12 

0.31 

0.53 

A75s 

0.75 

10-2 

9.3 • 1025 

1.5 

49.5 

1.5 

505 _ 1605 

0.49 

18 


1.0 

A75oi 

10-1 

9.3 • 1026 

0, 1.5 

1005 

0.16 

9 


0.29 

BIO 

1.0 







0.021 

0.11 

11 

18 

0.86 

1.5 

B75 

0.75 







0.0059 

0.031 

6.5 

11 

0.58 

1.0 

B50 

0.5 

10-2 

9.3 • 1025 

0.24 

7.8" 

0, 0.5, 1.5 

805 

0.0024 

0.013 

4.5 

7.9 

0.45 

0.78 

B25 

0.25 







0.0011 

0.0063 

3.4 

5.9 

0.36 

0.63 

BOO 

0.0 







0.00068 

0.0035 

2.7 

4.7 

0.31 

0.53 

B75^ 

0.75 

10-2 

9.3 • 1025 

0.24 

7.8 

0.5-1.5 

805 

0.0059 

-0.031 

6.5- 

11 

0.58-1.0 

B75oi 

10-* 

9.3 • 1026 

0, 1.5 

0.0097 

3.6 

0.32 

CIO 

1.0 







0.0043 

0.022 

6.3 

11 

1.5 

0.86 

C75 

0.75 







0.0012 

0.0063 

3.8 

6.6 

1.0 

0.58 

C50 

0.5 

10-2 

9.3 • 1025 

0.082 

2J‘‘ 

0, 0.5, 1.5 

805 

0.00049 

0.0026 

2.7 

4.6 

0.78 

0.45 

C25 

0.25 







0.00025 

0.0013 

2.0 

3.5 

0.63 

0.36 

COO 

0.0 







0.00014 

0.00072 

1.6 

2.8 

0.53 

0.31 

SSMol 


10-‘ 

9.3 • 1026 



1.05 


0.012" 

3.9 

0.3 

SSM 02 

0.84 

10-2 

9.3 • 1025 

0.3 

10 

0.48 

805 

0.012" 

8.5 

0.65 

SSE 04 


10-4 

9.3 • 1025 



1.98 


1.0 

171 


13 


Notes. All simulations are initially perturbed by a small plane sinusoidal wave in the v-direction as in FP2015. 

Considering of = H 2 and/5 = He. 0 = Earth T = 1.5rc,He. T ^ Tcmb- M = Mmooh- ^ 7j = 0.5. yj = 1.5. 



Fig. 6: Two-dimensional density map of three-dimensional 
space. 


by the Lennard-Jones potential and planetoids for clumps that 
are principally bound by gravity. If a fluid is in a phase transi¬ 
tion, it is important to distinguish between a strong gravitational 
potential above the ideal gas Jeans criterion and a weak grav- 



Fig. 7: Colour mapping of density map. 


itational potential below it. In the strong gravity case, a grav¬ 
itational collapse happens, leading to the formation of a hot, 
gaseous planetoid. Phase transitions only happen at the begin¬ 
ning as the fluid heats up above the critical temperature where 
no solid comets can form. 

In the weak gravity case, no gravitational collapse happens 
and solid comets form thanks to the phase transition. The comets 
attract each other gravitationally, which leads to the formation 
of a solid planetoid. During the comet aggregation, the number 
of bound molecules does not rise. This means that the plane¬ 
toid only captures comets and no single molecules. Therefore, 
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the fraction of bound molecules is identical to the number of 
molecules that underwent phase transition. 

In this article, we compare fluids with different molecu¬ 
lar fractions by keeping the physical properties alike (con¬ 
stant TITc,a and nlric^a)- By decreasing Xq, the mean mass per 
molecule m increases, since rria < nip. It is therefore not possi¬ 
ble to have the same number of molecules per super-molecule rj, 
the same super-molecule mass Msu = nip, and the same gravity 
Gj oc at the same time. 

Since we want to study the reaction of fluids with different Xq, 
above and below the ideal gas Jeans criterion, we need to ensure 
that 7 j is > 1 or < 1 in the compared simulations. For that reason, 
neither the total mass nor the number of molecules remains the 
same when comparing fluids with different Xq,, while yj remains 
the same. Since G 2 (Equ. and p^^^ oc nT^G (Equ. 39) , 
then p 2 nT^ and thus M = mp 2 nT^ . Therefore, changing m 
by a factor 2 leads to a total mass factor decrease of 32. 

Different cases are studied: T = T = l.STc^p and 

T = (7’cMB/7’c,H2)7’c,a with n = 10"^nc,a and = 1, 0.75, 
0.5, 0.25, and 0. These initial parameters are shown in Eig. 3. 
In addition, different densities are used for simulations B75 and 
C75 to compare cases with n > ri- with n < ri- (see Sect. 2.3.3) . 
We also study solar system abundances for different total masses 
and number densities. The simulations are summarized in Table 
i. 

The simulations with yj < 1 were run until a steady-state so¬ 
lution was reached. Most simulations with yj > 1 were stopped 
once a planetoid formed, which typically happens after ~5t. 
Even if the resulting fluid did not reached a steady state then, 
no further developments are expected, as checked in EP2015. A 
few simulations were run for 15t to confirm this. At ^ = 5t, 
the planetoid of B75 with y = 1.5 is 15% hotter than average, 
whereas unbound molecules are 5% colder. After another lOr, 
however, at ^ = 15t, the temperature differences are below 1%. 
Similarly, at r = 5 t, there are very strong regional temperature 
differences of > 10% considering subdomains of 1% volume. At 
t = 15t, they are < 5%. 

4.1. Above critical temperatures 

In this Section, we consider fluids where both Tq. and Tp are 
above the critical temperature. The number densities correspond¬ 
ing to the different abundance number fractions can be seen in 
Fig. 3. 

4.1.1 . Time evolution 

Pigure_8 shows the global temperature evolution of the A simula¬ 
tions (see Table 1). In all simulations, the weakly self-gravitating 
fluids with yj = 0.5 and the fluids without gravity are very simi¬ 
lar and do not react significantly to the velocity perturbation. On 
the other hand, as expected for the sufficiently self-gravitating 
fluids with yj = 1.5 the introduced perturbation rises exponen¬ 
tially. The reaction time is similar for all number fractions, but 
the mixtures reach slightly higher temperatures. Keeping in mind 
that all simulations have the same yj value, but the total mass dif¬ 
fers with Ml I M 2 = (jnilni 2 )~^. 

To get a deeper understanding of the internal processes, we 
need to distinguish between the a and p component. Eigure ^A 
shows the fraction of bound molecules in the simulations as a 
function of time. No phase transition happens above the critical 
temperature, therefore no comets form in the simulations with 
yj = 0 and 0.5. With yj = 1.5, the gravitational collapse leads to 



Eig. 8: Global temperature as a function of time of the A simu¬ 
lations with yj = 1.5 and 0.5. 



0 10 20 30 40 50 


R K] 

Eig. 11: Density of planetoid as a function of the radius of the 
simulations A75 and B75 with yj = 1.5at^ = 5T and B75 with 
yj = 0.5 at r = 30t. 


the formation of a planetoid. In its centre, the gravitational pull 
is strong enough to keep the molecules bound even though the 
temperature is well above the critical value. 

The yS-molecules, as they are twice as heavy as the a- 
molecules, fall faster into the forming planetoid. Indeed, even 
in A75, with only xp = 0.25, the fraction of bound yS-molecules 
surpasses the fraction of bound a-molecules for r > 2. 

4.1.2. Planetoid formation 

Eigure 10a (page 20) shows a time sequence of snapshots and 
super-molecule, comet-size distributions condensed as grains or 
comets. The parameter Acomet is the number of super-molecules 
in one comet and /(Ab) is the comet size distribution function. 
At the beginning with t < 3t, small comets of either a- or p- 
molecules form. Ait = 3t, a planetoid with A O.lAtot forms 
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7j = 1.5 7j=0.5 



Fig. 9: Evolution of the number of bound molecules as a function of time of the A, B, and C simulations with yj = 1.5 and 0.5. 
The simulations with yj = 0.5 are stopped after 5t, 30t, and 20t for the A, B, and C simulations, respectively, when no further 
significant development is expected. 


consisting of both components. One can already see a dominance 
of yS-molecules, especially in the centre. 

Beginning at ^ = 3t, and even more clearly at ^ = 4t, one 
can observe the formation of a big core consisting only of (3- 
molecules (isolated yS-dot). In the snapshots this corresponds to 
the planetoid shown as a /5-core surrounded by cr-molecules. 
Once the planetoid has reached this form, it reaches a steady 
state. Its temperature matches the gas temperature, and the tem¬ 
perature fiuctuations level out (see FP2015 for more details on 
planetoid and comet temperatures). 

Figure JJ_ (top) shows the planetoid density of simulation 
A75 as a function of radius. Even though the fiuid consists 
of only 25% yS-molecules, the planetoid consists mostly of (3- 
molecules with fp = 0.86. The a-molecules are only a small 
fraction and mostly present in the outer part. The gaseous nature 
of this body is visible as the density regularly decreases in radii. 

4.1.3. Scaling 

The scaling of simulations using super-molecules has already 
been discussed in FP2015. In order to obtain the correct be¬ 
haviour, the gravitational forces Fq need to be small on inter- 



Fig. 12: Fraction of bound molecules as a function of time for 
the simulations A75s with yj = 1.5 and different A/tot • 


molecular scales compared to the Fennard-Jones forces Flj, i.e. 
^g(^c) < ^lj(^c )5 where is the cut-off radius (see Equ. 40). A 
turning point can be identified up to which (A^comet) follows 
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Fig. 13: Fraction of boundyS-molecules as a function of A^tot for 
the simulations A75s. 


the power-law Nb (A^comet) ^ with negative index < -1, 

whereas after the turning point Nb (A^comet) follows a second 
power law with index = I (see Fig. 10a). The appearance 
time of this turning point is independent of Atot, whereas the size 
of the comet at the turning point scales as Acomet/^^tot ~ 
thus Acomet ~ 10. This corresponds roughly to the smallest num¬ 
ber of nearest neighbours in the condensed phase in 3D for which 
surface effects start to be dominated by volume effects. 

Figure 12 shows the fraction of bound molecules as a func¬ 
tion of time for all A75s simulations with Atot = 50^ to 160^. As 
shown in FP2015, the slight time delay between the simulations 
can be attributed to the random seed. In any case, the asymptotic 
final value is physically more important, and is the same for all 
Atot when considering both components. The final value of the 
yS-molecules, on the other hand, very slightly declines with in¬ 
creasing Atot as can be seen in Fig. JA. It follows the power-law 
NB/Ntot ~ 0.213 over the range Atot = 10^ - 10^ ^ 

4.1.4. Extrapolation to physical scale 

The simulations should actually represent a H 2 -He fluid mixture 
with ~10^^ molecules. As outrageous as this extrapolation might 
appear, this is exactly what usually takes place in many other 
types of simulations (cosmological, galactic, or stellar simula¬ 
tions) because as long as the physical scale invariant aspect of 
the physics between the macro- and micro-scales are separated 
by enough orders of magnitude the exact range of scale differ¬ 
ence does not matter over dynamical timescales. For longer sim¬ 
ulation timescales one can check how the results scale with A by 
running simulations with different A, which is the reason why 
we always run the simulations with several A. Extrapolating the 
previous power law to physical scales, we find that ~3% of (3- 
molecules settle inside the planetoid, instead of -15%. Thus the 
simulations overestimate species segregation, which is to be ex¬ 
pected in view of the increased fluctuations when the number 
of particles decreases. Segregation effects should be treated with 
caution, as we are extrapolating values in a range that is less 
than two orders of magnitude or 45 orders of magnitude away. 
Larger simulations should allow us to better constrain the effec¬ 
tive species segregation in realistic conditions. 

4.2. Between critical temperatures 

In this section, we consider fluids with and Tp > 

with different component fraction Xa. The number density n has 
been chosen in such a way that for the molecular fractions Vq, > 
0, the Qf-component with number density = Xq, • ^ is in a phase 
transition. 


As the temperature of the B simulations is an order of mag¬ 
nitude smaller than in the A simulations, the same is the case 
for the gravitational potential (see Equ. ^). Eor that reason, it is 
sufficient to use Atot = 80^ for these simulations. 

4.2.1. Above the ideal gas Jeans criterion 

Eigure 9B on the left side shows the time evolution of the fraction 
of bound molecules for the B simulations with yj = 1.5. One 
can see the similarity to Eig. 9A, but the fluids with a high Xq, 
value are rising to higher values even before the perturbation is 
becoming dominant. This is because the a-portion of the fluid 
is in a phase transition and small ice grains are forming even 
without the help of gravity. The formed planetoid is gaseous, 
as can be seen in Eig. jj_ (middle). This shows that the phase 
transition does not have an important effect if yj > 1 and that the 
instability can be predicted by the ideal gas Jeans criterion. 

The density at the core of the planetoid of simulation B75 
is lower than that of A75. This is explained by the fact that by 
keeping yj = 1.5, the value for Gj is lower for the B simula¬ 
tions than for the A simulations as Gj oc T (see Equ. 45). Having 
a lower gravitational potential, the density at which the repul¬ 
sive Lennard-Jones term and the attractive gravitational term are 
equal is lower. 

The planetoid is still dominated by /J-molecules, but there are 
more tr-molecules than in the A75 simulation with fa = 0.26 and 
fp = 0.74. As in the A simulations, the core consists of gaseous 
He, as is clearly visible in Eig. 10b (page 20). 

4.2.2. Below the ideal gas Jeans criterion 

As can be seen in Eig. X the tr-component of the simulations 
BIO, B75, B50, and B25 all lie on the Maxwell line and are thus 
in a phase transition, which implies, according to Equ. (2), that 
they are gravitationally unstable even with yj < 1. 

The right side of Eigure 9B shows the evolution of the frac¬ 
tion of bound molecules of the B simulations with yj = 0.5. 
The timescale is much larger (30t instead of 5t in the case of 
yj = 1.5), having a smaller gravitational potential, the long- 
range gravitational term is lower, and therefore the creation of 
any potential comet or planetoid takes more time. 

The one-component fluids, consisting of either uniquely a- 
molecules (BIO) or/J-molecules (BOO) have already been studied 
in detail in EP2015. The a-fluid B10 is unstable as it is in a phase 
transition, whereas the yS-fiuid BOO is stable as its temperature is 
above the critical value and no phase transition is possible. 

The simulations of fluid mixtures B25, B50, B75 with yj = 
0.5 are all unstable, even gravitationally. We can distinguish a 
clear difference in the simulations with yj = 1.5 in that only 
the a-molecules form comets, whereas the /J-molecules remain 
in gaseous form. Even in the simulation B25, which has only 
25% Qf-molecules, the comets and planetoid consist almost ex¬ 
clusively of Qf-molecules. This difference between yj = 1.5 and 
0.5 can also be seen when comparing Eig. 14a with Eig. 10b 
(pages 21 and 20). 

Eigure JJ_ (bottom) shows the radius of the planetoid at 
t = 30t of B75 with yj = 0.5. Comparing with the planetoid 
of B75 with yj = 1.5, we see that the high-gravity planetoid 
consists mostly of /J-molecules in gas phase, whereas the low- 
gravity planetoid consists of mostly of-molecules in solid phase, 
surrounded by an atmosphere. Very few /J-molecules have been 
trapped during the planetoid formation, providing an interesting 
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Fig. 15: Fraction of bound yS-molecules as a function of time 
for the simulations B75y. The simulations are stopped once they 
reach asymptotic values. 


example of a body forming with a distinct composition from the 
original medium as a result of the initial phase transition state. 


4.2.3. Different yj values 

The previous sections show that a fluid in a phase transition 
above the ideal gas Jeans criterion, i.e. with yj = 1.5, forms 
a gaseous planetoid consisting mostly of yS-molecules due to a 
classical ideal gas Jeans collapse. On the other hand, a fluid in 
a phase transition with yj = 0.5 forms small cr-comets due to 
the phase transition. These comets are attracted to each other by 
gravity, leading to the formation of a rocky planetoid, consisting 
almost exclusively of cr-molecules. In this section, we vary yj 
from 0.5 to 1.5. 

Figure 15 shows the fraction of bound yS-molecules. It is ris¬ 
ing steeply for fluids with yj > 1 in accordance with the ideal gas 
Jeans criterion and the formed planetoid is gaseous and consists 
mostly of yS-molecules. The fluid with yj = 1 also produces a 
gaseous planetoid, but the percentage of /J-molecules is already 
dropping a little. Interestingly, in the fluids with 0.7 < yj < 1, 
the p fraction is also rising. The instability criterion of Equ. (2) 
is for all components, not only one of them. 

Figure 14b (page 2j_) shows snapshots and comet-size dis¬ 
tributions of simulation B15y with yj = 0.8. One sees that at 
first {t < 8t) only the a-molecules are collapsing and forming 
a rocky planetoid. Then, owing to the great attractive force of 
the Qf-planetoid, many yS-molecules gather around it, forming an 
atmosphere {t = 12t). A /J-atmosphere can also be observed, in 
a less striking way, for the simulation with yj = 0.5 in Fig. 14a. 
What happens afterwards is very interesting: at r = 16 t, one sees 
that the rocky planetoid swaps the a- and yS-molecules and the 
heavier yS-molecules replace the tr-molecules near the centre. 


4.3. Below critical temperatures 

In this Section, to complete the study of binary fluid mixtures, 
we consider fluids where both Ta and Tp are below the critical 
temperature. The number density n has been chosen in such a 
way that for the molecular fractions Vq, > 0, the a-component 
with number density = Xq, • ^ is in a phase transition. 



Fig. 17: Fraction of bound molecules with cluster mass md > 
mHe as a function of time of simulations in and out of the unviri- 
alizable density domain. 


4.3.1. Above the ideal gas Jeans criterion 

The left side of Figure (page shows the time evolution of 
bound molecules of the C simulations with yj = 1.5. There is a 
distinct difference compared to the A and B simulations, which 
form yj-planetoids; only the percentage of bound tr-molecules 
rises and the forming planetoid only consists of tr-molecules (see 
Fig. 16, page 22). This is slightly counter-intuitive at first, as one 
could expect the yS-molecules to be even more eager to fall into 
the planetoid than in the A and B simulations, since the temper¬ 
ature is lower. 

Owing to the very low temperature of the C-simulation, how¬ 
ever, the tr-molecules quickly form comets from the very begin¬ 
ning. These comets are heavier than the /J-molecules and decel¬ 
erate faster into the planetoid. 

4.3.2. Below the ideal Gas Jeans criterion 

The evolution of the simulations below the ideal gas Jeans cri¬ 
terion is analogous to the B simulations. The fraction of bound 
Qf-molecules in the pure a-fluid and the mixture rise, and the 
fraction of bound molecules of the pure yS-fluid remains very low. 
This is in accordance with Fig. 3 where the a-molecules are un¬ 
stable but theyj-molecules are stable. The simulations CIO, C75, 
C50, and C25 form a rocky a-planetoid, as already seen in the B 
simulations (see Fig. 14a) . 

4.4. Virial theorem 

When comparing the simulations above the ideal gas Jeans in¬ 
stability, there is a clear difference between the A and B simula¬ 
tions on one side, and the C simulations on the other. A gaseous 
yS-planetoid forms in the first two, whereas a rocky a-planetoid 
forms in the latter. Booking at the virial terms of the fluids (see 
Sec. 2.3.3) , Equ. (31) is fulfilled in the A and B simulations, 
whereas for the C simulation, the density is in the unvirializable 
domain D. In this Section, we vary the densities of the A and 
B simulations in order to be in and out of the unvirializable do¬ 
main. 

Figure shows the time evolution of clusters that have a 
higher mass than one /J-molecule (Nc\,a > 2) for the simulations 
in the unvirializable domain (A75oi and B75oi) and below (A75 
and B75). A very quick rise of H 2 comets for the unvirializable 
fluid happens, both with and without gravity, which is in ac- 
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Fig. 18: Density of planetoid at r = 5t as a function of the ra¬ 
dius of the simulations in and out of the unvirializable density 
domain. 



Fig. 19: Fraction of bound molecules as a function of time of 
the simulations B75, B75 with removedyS-molecules and BIO. 
Tj = 0.5. 


cordance with Equ. (31), as neither the repulsive Lennard-Jones 
term nor the kinetic energy can withhold the attractive Lennard- 
Jones term thus leading to the formation of comets. Even in sim¬ 
ulation A75oi, with a temperature above the critical temperature, 
this comet formation is taking place, even though a phase transi¬ 
tion is officially not possible. A slow comet formation only takes 
place for the virializable fluids. 

Once the exponential growth of the perturbation becomes 
important (t > 2t), the unvirializable fluids have created an im¬ 
portant number of comets heavier than the yS-molecules, which 
fall faster in the forming planetoid as a result of dynamical fric¬ 
tion. This can be seen in Eig. 18 where the planetoids of the sim¬ 
ulations A75 and B75 consist mostly of yS-molecules, whereas 
the planetoid of B75oi consists mostly of tr-molecules. A some¬ 
what special case is A75oi, where the planetoids composition is 
almost perfectly fffty-fffty. This can be explained by the fact that 
because is is above the critical temperature, the comets are not 
really solid, but consist of a dense gas that is able to mix easily 
withyS-molecules. Thus, once a a-planetoid has formed using all 
the heavy a-comets, the yS-comets fall into the planetoid and mix 
with it. 

4.5. Influence ofp-molecules on a-molecules below the ideal 
gas Jeans criterion 

As can be seen in Eig. 9C, almost no yS-molecules form comets 
if yj = 0.5 and the percentage of yS-molecules in the planetoid 
is negligible. Granted, the concentration of He around the plan- 



10 ® 10 '® 10 '’ 10 '® 10 '^ 10 '^ 10 ° 


Eig. 20: yj of different total ffuid masses at T = 10 K, as a func¬ 
tion of number density, indicating either ideal gas Jeans instabil¬ 
ity, or instability owing to phase transition. 

etoid rises slightly as can be barely seen in Eig. 14a. Thus the 
question can be raised whether a small fraction of a secondary 
molecule (such as He in the case of molecular clouds) needs to 
be included in low-gravity simulations. To answer that question, 
simulation B75 with yj = 0.5 was run again, but all yS-molecules 
were removed and their mass was equally distributed to the a- 
molecules to maintain the same gravitational potential. 

Eigures 19 shows the time evolution of the fraction of bound 
molecules of the simulations B75, B75 without yS-molecules and 
BIO for comparison. Even though the two B75 simulations are 
similar, there are differences to be observed. The fraction of 
bound Qf-molecules of the simulation B75 should correspond to 
the total fraction of bound molecules of the simulation without 
yS-molecules, but the latter is higher; the/J-molecules in B75 have 
a damping effect on the comet formation. In addition, the simu¬ 
lation without yS-molecules is rising to a higher value at the end 
of the simulation. 

The inclusion of a small fraction of a secondary molecule 
does change the look of the simulation by damping the comet 
formation of of-molecules. Eor that reason, the inclusion of sec¬ 
ondary molecules in more realistic simulations is useful. 

4.6. Physical systems 

Up to now, we have looked at theoretical models, varying Xq, 
from 0 to 1, and setting the temperature and density as a fraction 
of the respective a critical values. The critical values for H 2 are 
Tc = 32.97 K and = 9.34 • 10^^ m“^. In astrophysical condi¬ 
tions, the He mass fraction is between w;He,ss = 0.2741 for the 
solar system (Lodders 2003) and w;He,MW = 0.2486 for the ini¬ 
tial Big Bang mixture (Cyburt et al. 2008), which translates to 
number fractions XHe,ss = 0.1598 and XHe,MW = 0.1428. 

Eigure 20 shows yj as a function of the number density for 
solar system abundances (x = 0.16) and T = lOK with total 
masses equal to the Moon, Earth, Jupiter, and Sun. H 2 is then 
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Fig. 21: Snapshots and comet-size distributions of the simula¬ 
tions SSMOl, SSM02, and SSE04. The slice selects in depth 
20% of the super-molecules. The squares in the two lower left 
frames are the same size as the next upper frame. 


in a phase transition for ^ > 4 • 10^^ m“^; only a Moon mass or 
below can be in a phase transition and below the Jeans criterion. 
The fluid is unvirializable for n> 6 • 10^^ m“^. 

If we go to a lower temperature, say the CMB 2.7 K, a H 2 
phase transition takes place for n > 10^^ m“^. In that case, fluids 
with Earth mass would be chemically unstable below the ideal 
gas Jeans criterion for ^ < 3 • 10^^ m“^ and with Jupiter mass for 
n < 3 • 10^^ m“^. Eluids with Sun mass, on the other hand, cross 
7 j = 1 only in the gaseous phase of H 2 . The lowest unvirializable 
density = 6 • 10^^ m“^ does not change a lot with temperature. 

The number of EET mesh cells Afft ^ and the simulation 
timescale t oc L both directly depend on L cx 
total calculation duration scales as 4im oc r • A^fft oc Eor 

that reason, simulating a fluid at CMB temperature with densities 
below 10^^ m“^ would translate to extremely long simulation run 
times with today’s computers. In addition, the upper limit for 
the mass of super-molecules is msM,max ~ 5 • (see Equ. 

40). Thus, the minimum number of super-molecules A^tot,min = 
'^^SM,max IS ~2 • 10^, 6.5 • 10^, 6.5 • 10^^ for simulating an 
Earth, Jupiter, and Sun mass, respectively. Eor that reason, for 
the time being we content ourselves to studying systems up to 
total mass comparable to the Earth mass. 


4.6.1. Planetoid formation 

Three simulations were run at a temperature of T = 10 K, which 
is above the critical temperature of He and below that of H 2 , and 
thus in a similar regime as the B simulations. Two simulations 
have a total mass equal to the Moon, with n ^ 10^^ m“^ which is 
above the ideal gas Jeans criterion and in the unvirializable do¬ 
main, and n ^ 10^^ m“^, which is below the Jeans criterion, and 
one has a total mass equal to the Earth and with ^ 10^"^ m“^, 

which is above the criterion. The simulation parameters are given 

in Table X 

Eigure ^ shows the snapshot and comet-size distribution of 
the three simulations after the formation of a planetoid. The fluid 
of SSE04 is above the ideal gas Jeans criterion and we observe 
the formation a He-planetoid, surrounded by H 2 , similar to Sect. 

4.2.1. The evolution of simulation SSM02, which is below the 
ideal gas Jeans criterion, leads to the formation of a rocky H 2 
planetoid, similar to Sect. 4.2.2. 

In the case of SSMOl, the density lies in the unvirializable 
domain, resulting in a formation of many H 2 -grains that are 
heavier than the He-atoms from the very beginning. This leads 
to the formation of a H 2 -planetoid similar to Sect. 4.3.1. 

5. Conclusions 

In our first article, EP2015, we studied the gravitational instabil¬ 
ity of a fluid in a phase transition. We extrapolated the results to 
the ubiquitous H 2 and showed that the formation of cold, mostly 
undetectable comet- and even planet-sized rocky H 2 clumps is 
very possible. The use of only one component gives a good 
first impression, but in cosmic gases, there is a mass fraction 
of u; 25 + 2 % He atoms. 

In the present work, we studied binary fluid mixtures ana¬ 
lytically and via numerical simulations. The results show that, 
depending on the circumstances, either He or H 2 planetoids can 
form. 

5.1. Analytic results 

The stability of a multicomponent fluid mixture has already been 
studied in the literature, mostly to study fluid binaries consisting 
of baryonic and dark matter. The wave number below which a 
fluid mixture is unstable is the sum of the Jeans wave-numbers 
of each component. Since the Jeans wave number is inversely 
proportional to (dPldp)~^, which is equal to zero in the case of a 
phase transition, a fluid mixture is unstable as soon as one of 
its components is in a phase transition. Physically what hap¬ 
pens is that when one species is in a phase transition, an over¬ 
density only increases its condensed phase fraction at constant 
pressure, instead of increasing pressure and producing no global 
force to counter gravity. The transformation from the gas to the 
condensed phase continues until the species is fully condensed. 

We studied the evolution of unstable fluid mixtures with the 
widely used Lennard-Jones intermolecular potential, which re¬ 
produces the H 2 phase transition very well (but it reproduces the 
He transition, which is not essential in this work, less well). We 
showed, using the virial analysis of Lennard-Jones fluid mix¬ 
tures, that there is a unvirializable density-domain D within 
which the attractive forces dominate the repulsive forces for any 
total mass M and no virial equilibrium is possible. These states 
can be reached in strongly dynamical situations (e.g. during col¬ 
lapses) and are able to produce condensed comets particularly 
quickly. Dynamical friction is important to separate species and 
condensed comets. Eor instance, if H 2 is in a phase transition, the 


Article number, page 13 of 22 



















A&A proofs: manuscript no. FP2016 


formed H 2 comets are heavier than the He-molecules, and pre¬ 
cipitate in a gravitational field, producing almost pure H 2 bodies. 

There are three reasons to concentrate on plane-parallel ini¬ 
tial collapses, as described in more detail in App. B: 

1. In typical cosmic conditions, the fastest collapsing geometry 
is sheet-like, not filament- or point-like. 

2. The adiabatic matter compression during collapse leads to 
the least heating in sheet-like geometry: in a sheet-like adia¬ 
batic collapse the gravitational energy released to the fluid is 
finite and amounts to a maximum increase of temperature by 
only a factor of about two, while in filament-like collapses 
the temperature diverges logarithmically as a function of fil¬ 
ament radius, and in point-like collapses the temperature di¬ 
verges as the inverse sphere radius. 

3. Radiative cooling is the easiest in sheet-like collapse. Indeed 
the absorption probability in sheet-like geometry remains al¬ 
most unchanged for any compression, and an initially trans¬ 
parent medium remains transparent, whereas the probability 
converges to one in filament-like and point-like geometries. 
Therefore, radiative cooling is barely slowed down in sheet¬ 
like collapses and, unlike in spherical or filament collapses, 
opacity is unable to prevent density from reaching high val¬ 
ues. This is a crucial point for this study, as the ISM con¬ 
ditions are commonly thought to be far away from the H 2 
phase transition conditions. 

5.2. Simulations 

As in FP2015, we used super-molecules to combine the Lennard- 
Jones intermolecular potential together with the gravitational 
potential in numerical simulations. Several binary fluid mix¬ 
tures were studied using two components: a and fi. Their re¬ 
spective properties (the most important being mal^np < 1 and 
Tc,alTc,p >1) were chosen to mimic a H 2 -He fluid, but the gen¬ 
eral properties of the fluids were made molecule independent. 

Three temperature domains can be defined: (A) above both 
critical temperatures, (B) between the critical temperatures, and 
(C) below both critical temperatures. In all three cases, the 
molecular fraction was varied and the fluids were simulated 
above and below the Jeans criterion. We used different numbers 
of molecules to test the scaling of the simulations. The precise 
number of super-molecules is not important for dynamical pro¬ 
cesses, but we found a weak dependence for segregation effects 
in the sense that coarser simulations exaggerate these effects. 

In case (A), both components are gaseous and an introduced 
perturbation does not grow when the gravitational potential is 
below the Jeans criterion. When above the Jeans criterion, the 
fluid collapses and forms a gaseous planetoid. The yS-molecules 
are twice as massive as the a-molecules, and fall faster into the 
planetoid. For that reason, the planetoid consists mostly of fi- 
molecules, surrounded by an a-atmosphere. This is independent 
of the molecular fraction even at very high Xq,- values, the 
planetoids consists mostly of yS-molecules. 

In case (B) and (C), the number density of the fluids was 
chosen so that the a-component is in a phase transition for all 
Xq, > 0. In fact, both cases are very similar since in both cases 
the yS-component is not in a phase transition. When the fluids are 
below the Jeans criterion, an instability happens because of the 
phase transition of the a-component, which leads to the forma¬ 
tion of H 2 comets and ultimately a rocky a-planetoid. This plan¬ 
etoid is surrounded by a yS-atmosphere, which is getting more 
important with increasing gravitational potential. As in case (A), 
the molecular fraction Xq, does not matter, even at very low 
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Fig. 22: Gravitational instabilities at different temperatures and 
densities for a fluid with Jupiter mass. 


XQ,-values, the planetoid still consists almost exclusively of a- 
molecules. 

A suprising observation occurs for cases (B) or (C) above 
the Jeans criterion. In that case, there is a race between the for¬ 
mation of small Qf-grains owing to the phase transition and the 
exponential growth of the perturbation. The heaviest bodies are 
decelerated faster and fall into the forming planetoid first. When 
the Qf-component is either gaseous or only forming very few 
and small comets, a yS-planetoid forms. On the other hand, if 
the Qf-component forms many grains that are heavier than the 
yS-molecules, an a-planetoid forms. We showed in the simula¬ 
tions that this race between a and fi is linked with the unviri- 
alizable density domain D. If a fluid reaches this domain, the 
a-component wins, otherwise the yS-component wins. 

5.2.1. Solar system abundances 

In addition to the above-mentioned simulations, fluids with so¬ 
lar system abundances and Moon or Earth mass were simulated. 
As shown in Fig. 20, a fluid with Earth mass cannot be below 
the Jeans criterion and still in a phase transition, but with Moon 
mass, this is possible. In that case, a rocky H 2 planetoid results. 
With a mass as low as the Moon, the fluid needs to be very dense 
to be above the Jeans criterion. In fact, the fluid would lie in 
the unvirializable density-domain D and, thereby, a H 2 -planetoid 
forms. For a fluid with Earth mass, on the other hand, even a rel¬ 
atively low-density fluid is still above the Jeans criterion. The 
result is a gaseous He planetoid with a H 2 atmosphere. 

5.3. Instability in H 2 -He fluid 

Figure 22 shows different possible planetoid and comet forma¬ 
tions due to gravitational instability for a fluid with Jupiter mass. 
A fluid is gaseous if it is below the phase transition domain and a 
fluid is solid or liquid if above. When the density is in the phase 
transition, it can rise without an increase of pressure. 

There can be no formation below the Jeans criterion if the 
fluid is not in a phase transition. Most of the planetoids due to 
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an ideal gas Jeans collapse consist of gaseous He, but if the fluid 
is in the unvirializable domain D, then a H 2 planetoid forms. 
This H 2 planetoid can be solid/liquid or gaseous depending on 
its temperature. If gaseous, He is able to percolate down, slowly 
transforming it into a He planetoid. 

If the fluid is in a phase transition, we have to distinguish be¬ 
tween a collapse above the ideal Jeans criterion, which leads to a 
gaseous He planetoid except in the unvirializable domain, where 
it becomes a rocky H 2 planetoid, and in a collapse below the 
ideal Jeans criterion, which also leads to a rocky H 2 planetoid. 

The usual average density domain of molecular clouds lies 
between 10^ and 10^^ H 2 /m^ and, with such a density, a H 2 phase 
transition is only possible at temperatures below ~5 K. However, 
molecular clouds are observed to follow a fractal mass distri¬ 
bution over a minimum of 4-6 orders of magnitude in column 
densities, so the average density is not a quantity to characterize 
molecular clouds properly. Since we know that stars form with 
densities H 2 /m^, by continuing this argument, intermedi¬ 
ate states covering all this density interval have to exist. 

Fluids with a high total mass, especially with stellar mass or 
above, reach the ideal gas Jeans criterion very quickly leading to 
gaseous He-planetoids. Fluids with lower total mass, however, 
as for example the cold globules observed in the Helix nebula, 
especially with Earth mass and below, have the ideal gas Jeans 
criterion at much higher densities and are in the phase transition 
domain before being above the ideal gas Jeans criterion. 

5.4. Perspectives 

This and the previous FP2015 study show that the cold ISM 
physics is much richer than previously imagined. The formation 
of substellar gaseous or rocky condensed bodies by the H 2 phase 
transition combined to gravity, appears natural once we recog¬ 
nize that collapses proceed most of the time along the sequence 
pancake, filament, and point, and in the first sheet-like phase 
high densities allowing a H 2 phase transition can be reached if 
the initial medium temperature is below ~15K. This tempera¬ 
ture limit would be even higher if radiative cooling had been 
considered. In the isothermal case this limit is ~33 K. 

Most of the ISM cold gas must therefore pass over molecular 
cloud lifetimes (~10^ -10^ yr) through such brief (~10^ -10"^ yr) 
singular sheet-like collapses where density diverges but not tem¬ 
perature. Observationally, such events are difficult to detect be¬ 
cause of the limited increase of temperature, opacity, and column 
density all along the collapse, while reaching high volume den¬ 
sities. When seen edge-on such sheet-like collapses would look 
like filaments. 

The simulations we were able to perform are still very lim¬ 
ited in total mass. Including He is necessary but this provides a 
number of complications with respect to the pure H 2 case, and 
widens the general picture found in FP2015. Combining the ac¬ 
cumulated experience of large-scale gas phase simulations by 
other authors (e.g. Renaud et al. 2013; Butler et al. 2015), we 
can easily extrapolate what larger simulations should produce 
with micro-AU resolution. Instead of one planetoid per simu¬ 
lation box, pc-sized sheet-like collapses should show filaments 
with longer lifetimes, which would funnel H 2 condensed bodies 
and produce a spectrum of planetoids, comets, and occasionally 
stars. The leftover condensed cold substellar bodies should then 
start to evaporate according to the ambient radiation flux and 
depth of their gravitational potential. The lifetime of such bod¬ 
ies should be short near the centre of galaxies, but much longer 
at the periphery of galaxies, or even in intergalactic space, es¬ 
pecially in cosmic filaments. One can postulate that, especially 


at the periphery of disk galaxies where the radiation heating is 
low, some fraction of the dark baryons can be trapped in the form 
of such condensed bodies. We plan to pursue further simulation 
work to deepen our understanding of the processes associating 
phase transition with gravitational dynamics. 
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Appendix A: Jeans instability 

We first recall the classical Jeans criterion for a one-component 
fluid, and then we show how the same approach can be used 
to find the solution of a two-component fiuid. See Grishchuk & 
Zeldovich (1981) for the solution of an /t-component fluid. 


Appendix A.1: One component 

The equations for conservation of mass and momentum and for 
the gravitational potential of a fiuid are written as 


dpv 

dt 


Following Jeans (1902), we supersede these equations with per¬ 
turbation terms in the v direction p = po Sp, P = Pq SP, 
O = Oo + dO and v = Vq Sv with Sv = (6v, 0,0)^, linearizing 
the equations and setting SP = 6p, i.e. 


f+V.(p.) = 

0, 

(A.1) 

+ V • (pvv) + VP + pVO = 

0 , 

(A.2) 

V • VO - 47rGp = 

0 . 

(A.3) 


d dp 

^-rpoV-di; = 0, 

dt 


dSv ldP\ 

Po-^+l^^j Vdp+poVdO 


0 , 


V • VdO - AnGdp = 0 . 


- iiop + ikpov 

= 0 , 

(A.7) 

— p + ikpQ<^ 

dp Is 

= 0 , 

(A. 8 ) 

— 4nGp 

= 0 , 

(A.9) 


which can be written in matrix form A • jc = 0, 


oj kpo 0 

k(dPldp)s cxJpo kpo 

-4nG 0 -k^ 



P 


0 


V 


0 


6 


0 


Non-trivial solutions for x require that the determinant of A van¬ 
ishes. 


det(A) = k^po 


dP 

(j)^ -h 4nGpo 


OJ' 




dp 


which is the case for either k = 0 , or 
AnGpo • 

A fiuid is unstable if o/ < 0, which is the case if 
< k] = 


72 ^ r2 _ 4;rGpo 


Appendix A.2: Two components 

Having two components a and p, the mass and momentum con¬ 
servation have to be fulfilled for each component as follows: 


^+^-(PaVa) = 

dt 

0 , 

(A. 14) 

^PB 

-^+^-(ppVp) = 

0 , 

(A.15) 

+ V • (Pq,i;q,i;q,) + VPq,+ pVO - 
dt 

0 , 

(A.16) 

+ V • (ppVpVi^ + VPp + pVO — 

0 , 

(A.17) 

V-V<D-47rG(p„+p/3) = 

0 . 

(A. 18) 


Superseding, as in App. A.l, these equations with perturbation 
terms A^ = A^o + dA^ and Ap = Apo + ^Ap in the v direction and 
linearizing them yields 


ddpa 

dt 

ddpp 

dt 


+ Pa^ ’dVa 
+ ppV • dvp 


(A.4) 

(A.5) 

(A. 6 ) 


ddVa IdPA 

Par—VdpQ,+PQ,VdO 


ddv, 


P/5 




dP^ 




0 , 

0 , 

0, 

0, 

0, 


(A.19) 

(A.20) 

(A.21) 

(A.22) 

(A.23) 


This system of partial differential equations is transformed to 
an algebraic system of linear equations in the Fourier space: 
dA = f dkAik) exp[/(^v - cot)], where A represents p, v, and O. 
The passage to Fourier space transforms the differential opera¬ 
tors djdt and djdx to multiplications by -ico and ik, respectively. 


V • VdO - 47rG(dpQ, + dpp) 
which transform into a linear equation system in Fourier space. 


— ioj pci + ik paO Vci 

= 0 , 

(A.24) 

-iojp/s + ikppovp 

= 0 , 

(A.25) 

-ioJpaOVa + ik I ^1 Pa + ikpaO^ 
\dpals 

= 0 , 

(A.26) 

IdPA 

-lOjppoVp + Ik j Pjs + ikppo^ 

= 0 , 

(A.27) 

-k^^ - AnGipa + pf}) 

= 0 . 

(A.28) 


(A. 10) 


OJ 

0 

kpaO 

0 

0 


Pa 


roi 

0 

OJ 

0 

kppo 

0 


h 


0 


0 

^PaO 

0 

k PaO 


Va 


0 

0 

kcl 

0 

OJppo 

kppo 


h 


0 

-47rG 

-AnG 

0 

0 

-k\ 


6 


0 


{dPldp)s 


(A.11) 


(A. 12) 


(A.13) 


This can be written in the matrix form A • jc = 0, defining 
{dPaldpa)s and c| = idPjjldp/ 3 )s as follows: 


(A.29) 


In order to simplify, we set = AnGpao and Tp = AnGppo and 
find the following determinant: 

det(A) = k^PaoPpo + + k^icl + cj)) a? 

+k^{-Vpcl-T„cl + k^clcfj\. (A.30) 

Again, to have a non-trivial solution, its determinant must be 
zero, which, in the case of k 9 ^ 0 , is 

+ (r„ + - k\cl+cl)) a? - e {Tpcl + - k^clcj) = 0, 

(A.31) 
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with the following solution: 



a fluid is unstable for < 0 or e 3, which is the case for 
^ ^Gz- 



Appendix A.2.1: Phase transition 

In the case of a phase transition, one of the pressure derivatives 
is equal to zero. Setting Cq, = 0 in Equ. (A.31) we get 

(j' + [r„ + u? - Tak^cj = 0 , (A.35) 

and its solutions is written as 



(A.36) 


Setting 0 ? = Equ. (A.35) , only the trivial ^ = 0 is a solution. 
Since 

(r« + < (r<, + + 4 Tak^cj, (A.37) 


Eig. B.l: Collapsing geometries. 


Appendix B.1: Gravitational energy 

The gravitational energy difference between the initial sphere 
and subsequent ellipsoids must be released as additional thermal 
energy. The gravitational energy of a revolution ellipsoid, with 
Eg,o = E’G,sphere(^) = -(3/5)GM^/r (Landau & Lifshitz 1975), 


is written as 


^oblate ('^) _ 

arccos^Z 

Eg,o 

Vl - Z-2 

-^prolate ('^) 

arcosh ( Vz) 

Eg,o 

Vl -z-i 

-^sphere ('^) 

z'^l 

Eg,o 



- - Z“^ + 0(Z“b, 

2 

(B.l) 

log(2Vz) + o|i^| 

, (B.2) 


(B.3) 


the upper sign solution of Equ. (A.36) is always positive and the 
lower sign solution is always negative for any k. Therefore one 
tu-solution of Equ. (A.36) is always negative and thus unstable, 
independent of the strength of either or 


Appendix B: Energy and radiation transfer during 
the contraction of a sphere towards an eiiipsoid 

We consider a non-rotating sphere of radius r initially in unsta¬ 
ble equilibrium, which contracts at constant mass as an ellipsoid 
with semi-principal axes a, b, and c (see Eig. B.l) . In a sheet¬ 
like collapse, two semi-axes remain the same (a = h = r) while 
one is decreasing (c = sr), leading to an oblate spheroid. In a 
filament-like collapse, one semi-axis remains the same (a = r), 
while two are decreasing together (h = c = sr), leading to a 
prolate spheroid. In a point-like collapse, all the three semi-axes 
decrease together (a = b = c = sr), remaining a sphere. Dur¬ 
ing compression, density increases by a factor Z = n/no. Since 
the ellipsoid volume is Ven = 4/3 nabc, compression changes as: 

^oblate ~ ^ 9 ^prolate — Z ^ , and ^sphere — Z ^ . 


Sheet-like contraction leads to infinite densities with finite tem¬ 
perature increase, which is much more favourable for reaching 
condensation conditions that filament-like or point-like contrac¬ 
tions. 

We show now that the maximum relative temperature in¬ 
crease of an initial perfect gas sphere initially in equilibrium is 
bounded. State 0 is the initial (unstable) equilibrium sphere case, 
and state 1 is any later, denser case that is not necessarily in equi¬ 
librium. Since in equilibrium, the initial state respects the virial 
condition, 

Eg,o + 2 £’kin,o = 0 9 (B.4) 

where £’kin,o is the kinetic energy. Since at rest, the initial sphere 
kinetic energy consists only of microscopic motion, and is pro¬ 
portional to the initial temperature Tq. 

The initial and later total energies are, 

-^tobO ~ ^G ,0 F ^kin ,0 9 (B.5) 

Etot,i = Eg,I + E’kin,! • (B. 6 ) 
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Taking into account possible radiative cooling, we suppose 
Etot,i ^ ^tot,o, which leads to, using the initial virial condition. 


Tl ^ £^kin,l ^ 2 ^G,l ^ 

To T’kin,0 ^0,0 


(B.7) 


The first inequality takes into account that state 1 is not neces¬ 
sarily in equilibrium; some kinetic energy may be attributed to 
macroscopic motion. 

Thus, using the above potential energy ratios, in the case of 
an oblate spheroid contraction, 

— <7T- I - 2Z“‘ + 0(Z“^), (B.8) 

To 


that is, the final temperature cannot exceed tt - 1 ^2.1 times the 
initial temperature. In the case of a prolate spheroid contraction, 
temperature is logarithmically bounded as Z increases. 


^ < log(4Z) -1+0 
^0 


log(Z) \ 

Z /’ 


(B.9) 


while in a spherical contraction, temperature is bounded by the 
cubic root of compression, 

T 

(B.IO) 

To 


Appendix B.2: Radiative cooling 

Energy lost by radiation lowers temperature, but if opacity in¬ 
creases during contraction at some point the radiative cooling 
rate drops below the heating rate as a result of gravitational en¬ 
ergy conversion, thereby slowing down the collapse. Here we 
show with simple arguments how opacity changes when con¬ 
tinuously contracting an initial sphere towards denser, smaller 
spheres, or towards denser revolution of oblate or prolate kinds 
of ellipsoids, keeping the longest axes constant and assuming 
uniform densities at any stage and constant absorption cross sec¬ 
tions. 


Appendix B.2.1: Optical depth 

The optical depth r in the cumulated absorption over a pho¬ 
ton path f: T = crndf , where cr is the absorption cross 
section of individual atoms with number density n. The cen¬ 
tral optical depth, calculated from the centre to the ellipsoid 
edge along some direction, is a first order estimator of the av¬ 
erage optical depth. We compare the optical depth tq = rcrriQ 
for the initial sphere with the later spheres. For revolution el¬ 
lipsoids, where a and c are the semi-long and short axes, re¬ 
spectively, the distance from the centre to some point on the 

edge is f(6) = ac/ V^2”sh?^T~c^~cos^ for oblate spheroids and 

i{6) = acj V^2~cos^~^T~c^”sin^ for prolate spheroids. The angle 
Q vanishes at the spheroid equator. Since the ellipticity s = cja 
varies as and Z“^/^ in the oblate, prolate, and spher¬ 

ical cases, respectively, the optical depth ratios as functions of 
compression Z and 6 are found to be 


T oblate 

1 

(B.ll) 

To 

Vsin^ 0 + Z~^ cos^ 0 

Tprolate 

Zl/2 

(B.12) 

To 

Vcos^ 0 + Z“* sin^ 0 

T sphere 

= z^i\ 

(B.13) 


To 



Fig. B.2: Absorption probability in contracting spheroids as 
function of compression Z > 1 calculated by Monte Carlo simu¬ 
lation. At Z = 1 all cases are spherical. The photons start either 
at the centre only or anywhere inside the ellipsoid in random di¬ 
rections. Different cases are represented where the initial sphere 
optical depth tq is indicated. 


Thus, in the oblate case the central optical depth ratio does not 
change along the poles and at high compression remains barely 
increased over most directions. In the prolate case it increases 
least along the equator, but is proportional to the square root of 
compression. In the spherical case it increases most rapidly as 
a power 2/3 of compression. Thus sheet-like compression pro¬ 
vides the least optical depth increase and spherical compression 
compression the most. 


Appendix B.2.2: Global absorption 

One can refine the previous estimate for cooling by calculating, 
for any point inside an ellipsoid, the probability for a photon to 
be absorbed. For a given optical depth r the absorption prob¬ 
ability is p = 1 - exp(-T). The global probability of absorp¬ 
tion must be calculated for all solid angles for all points. These 
4- or 5-dimensional integrals for bi- or tri-axial ellipsoids does 
not seem to be solvable analytically, and straightforward numer¬ 
ical quadratures would be expensive. Thus we resort to a Monte 
Carlo draw to estimate these quantities. We pick randomly and 
uniformly a number of points inside the ellipsoid and a random, 
uniform directional unit vector n, and find the two distances fi, 
^ 2 , to the edge of the ellipsoid, allowing us to calculate two op¬ 
tical depths Tl, T 2 , and the corresponding absorption probabili¬ 
ties pi, P 2 for each point. Knowing the starting position x inside 
the ellipsoid (a, b, c) and the direction vector n, we find the two 
signed distances to the ellipsoid edge by solving the quadratic 
equation (x inx)la^ + (?/ + + (z + jc^ = 1 for 

i. Explicitly, noting a = p = b~^, y = for each point 
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X = [x, y, z] the procedure becomes 

A = anl+Pnl + ynl , 

B = axrix + Pyriy + yznz, 

C = ax^ + /3y^ + yz^ - 1 , 


(B.14) 

(B.15) 

(B.16) 

(B.17) 

(B.18) 


D = V 52 - AC, 

= -{D + B)IA, €2 = (D-B)IA. 


For each set of average absorption probabilities can be found 
for a range of crs. The two absorption probabilities pi = I - 
Qxp(crn\^i\), i = 1,2, provide two distinct probabilities for each 
point. Each set of p/S should converge towards a similar average 
value. The difference allows us to check the error obtained with 
a finite number of points. Between 2-10"^ (sphere case) and 3-10^ 
points (oblate spheroid case) were drawn such that the log^o Pi 
between the two sets differ by at most 0.01. The result is shown 
in Fig. B.2. The error bars are comparable or smaller than the 
thickness of the curve. 

The sphere and prolate spheroid cases quickly become opti¬ 
cally thick, increasing as and respectively, in the op¬ 
tically thin regime. In contrast, the absorption of a contracting 
optically thin oblate spheroid increases logarithmically until it 
reaches Zto~ 1 beyond which it remains approximately constant; 
in other words if the initial state is optically thin, it remains so 
even after infinite compression. The emission signature of a col¬ 
lapsing sheet should therefore remain observationally barely no¬ 
ticeable, since both temperature and optical thickness increase 
by very modest factors in comparison with the other geometries. 

Fig. 2 shows how an initial sphere at T = 10 K, P = 10“^^ Pa 
would change its temperature and pressure when contracting adi- 
abatically, changing its gravitational energy into thermal energy. 
Clearly the sheet-like collapse is the most favourable geometry 
for reaching the H 2 phase transition regime. Including radiative 
cooling, which is the fastest in sheet-like geometry, can only help 
in this regard. 
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(a) A75 with yj = 1.5 (b) B75 with yj = 1.5 

Fig. 10: Time sequence of the simulations A75 and B75. On the left side, the slice shows in depth 20% of the super-molecules. On 
the right side, A^comet is the number of super-molecules in one comet and / (A^b) is the comet size distribution function. 
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(b) B75 with yj = 0.8 

Fig. 14: Time sequence of the simulation B75. On the left side, the slice shows in depth 20% of the super-molecules. On the right 
side, A^comet is the number of super-molecules in one comet and / (A^b) is the comet size distribution function. 
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Fig. 16: Time sequence of the simulation C75 with jq = 1.5. On the left side, the slice shows in depth 20% of the super-molecules. 
On the right side, A^comet is the number of super-molecules in one comet and / (A^b) is the comet size distribution function. 
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